Cardiac mechanics Benchmark - Problem 2

Cardiac mechanics Benchmark - Problem 2#

This code implements problem 2 in the Cardiac Mechanic Benchmark paper

Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.

from pathlib import Path
import dolfin
try:
    from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
    from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
from fenics_plotly import plot
here = Path(__file__).absolute().parent
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 here = Path(__file__).absolute().parent

NameError: name '__file__' is not defined
geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
# geometry = pulse.geometries.benchmark_ellipsoid_geometry()
2024-01-29 16:19:22,465 [925] INFO     pulse.geometry_utils: 
Load mesh from h5
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 10.0
material_parameters["bf"] = 1.0
material_parameters["bfs"] = 1.0
material_parameters["bt"] = 1.0
material = pulse.Guccione(parameters=material_parameters)
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(
        V,
        Constant((0.0, 0.0, 0.0)),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, lvp, 10.0, initial_number_of_steps=200)
2024-01-29 16:19:24,404 [925] INFO     pulse.iterate: Iterating to target control (f_28)...
2024-01-29 16:19:24,404 [925] INFO     pulse.iterate: Current control: f_28 = 0.000
2024-01-29 16:19:24,405 [925] INFO     pulse.iterate: Target: 10.000
2024-01-29 16:19:24,555 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,556 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,556 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,557 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:25,239 [925] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full
2024-01-29 16:19:25,449 [925] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full
2024-01-29 16:20:25,159 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,160 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,160 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,161 [925] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,319 [925] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full
2024-01-29 16:20:25,515 [925] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full
*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 44),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 82),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 107),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 140),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 173),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 255),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 296),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 337),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 378),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 419),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 460),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 501),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 542),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 583),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 608)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 42),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 80),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 105),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 138),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 171),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 253),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 294),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 335),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 376),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 417),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 458),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 499),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 540),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 581),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 606)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
endo_apex_marker = geometry.markers["APEX_ENDO"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
print(
    ("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
        endo_apex_pos[0],
    ),
)
Get longitudinal position of endocardial apex: 30.7991 mm
epi_apex_marker = geometry.markers["APEX_EPI"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
print(
    ("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
        epi_apex_pos[0],
    ),
)
Get longitudinal position of epicardial apex: 32.0417 mm
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()